Departamento de Física - Faculdade de Ciências e Tecnologia da Universidade de Coimbra

Física Computacional - Ficha 2 - Zeros de Funções

Rafael Isaque Santos - 2012144694 - Licenciatura em Física



In [1]:

    
import numpy as np
%matplotlib inline

Define-se a nossa função:

$f(x) = \sin (x)$

e a sua derivada:

$f'(x) = \cos (x)$



In [2]:

    
f = lambda x: np.sin(x)
df = lambda x: np.cos(x)
my_stop = 1.e-4
my_nitmax = 100000
my_cdif = 1.e-6

Método da Bisecção



In [3]:

    
def bi(a, b, fun, eps, nitmax):
    c = (a + b) / 2
    it = 1
    while np.abs(fun(c)) > eps and it < nitmax:
        if fun(a)*fun(c) < 0: b = c
        else: a = c
        c = (a + b) / 2
        it += 1
    return it, c, fun(c)



In [4]:

    
bi(2, 4, f, my_stop, my_nitmax)









    Out[4]:





(11, 3.1416015625, -8.9089102066436886e-06)

Método da Falsa Posição:



In [5]:

    
def regfalsi(a, b, fun, eps, nitmax):
    c = (a * fun(b) - b * fun(a)) / (fun(b) - fun(a))
    it = 1
    while np.abs(fun(c)) > eps and it < nitmax:
        if fun(a) * fun(c) < 0: b = c
        else: a = c
        c = (a * fun(b) - b * fun(a)) / (fun(b) - fun(a))
        it += 1
    return it, c, fun(c)



In [6]:

    
regfalsi(2, 4, f, my_stop, my_nitmax)









    Out[6]:





(3, 3.1415903579556947, 2.2956340984908623e-06)

Método de Newton-Raphson:



In [7]:

    
def newtraph(c0, fun, dfun, eps, nitmax):
    c = c0
    it = 1
    while np.abs(fun(c)) > eps and it < nitmax:
        c = c - fun(c) / dfun(c)
        it += 1
    return it, c, fun(c)



In [8]:

    
newtraph(2, f, df, my_stop, my_nitmax)









    Out[8]:





(6, 3.1415926536808043, -9.1011076112782062e-11)

Método da Secante



In [9]:

    
def secant(a, b, fun, eps, nitmax):
    c = (a * fun(b) - b * fun(a)) / (fun(b) - fun(a))
    it = 1
    while np.abs(fun(c)) > eps and it < nitmax:
        a = b
        b = c
        c = (a * fun(b) - b * fun(a)) / (fun(b) - fun(a))
        it += 1
    return it, c, fun(c)



In [10]:

    
secant(2, 4, f, my_stop, my_nitmax)









    Out[10]:





(3, 3.1415903579556947, 2.2956340984908623e-06)



In [11]:

    
def bi2(a, b, fun, eps, nitmax, cdif):
    c = (a + b) / 2
    c_prev = a
    it = 1
    while not((np.abs(fun(c)) < eps and np.abs(c - c_prev) < cdif) or it > nitmax):
        if fun(a)*fun(c) < 0: b = c
        else: a = c
        c_prev = c
        c = (a + b) / 2
        it += 1
    return it, c, fun(c), c_prev, fun(c_prev), (c - c_prev)



In [12]:

    
bi2(2, 4, f, my_stop, my_nitmax, my_cdif)









    Out[12]:





(21,
 3.1415929794311523,
 -3.2584135910528161e-07,
 3.141592025756836,
 6.2783295730092144e-07,
 9.5367431640625e-07)



In [13]:

    
def regfalsi2(a, b, fun, eps, nitmax, cdif):
    c = (a * fun(b) - b * fun(a)) / (fun(b) - fun(a))
    c_prev = c + cdif/2
    it = 1
    while not((np.abs(fun(c)) < eps and np.abs(c - c_prev) < cdif) or it > nitmax):
        if fun(a) * fun(c) < 0: b = c
        else: a = c
        c_prev = c
        c = (a * fun(b) - b * fun(a)) / (fun(b) - fun(a))
        it += 1
    return it, c, fun(c), c_prev, fun(c_prev), (c - c_prev)



In [14]:

    
regfalsi2(2, 4, f, my_stop, my_nitmax, my_cdif)









    Out[14]:





(5,
 3.1415926535897931,
 1.2246467991473532e-16,
 3.1415926536048882,
 -1.5094913867333564e-11,
 -1.5095036332013478e-11)



In [15]:

    
def newtraph2(c0, fun, dfun, eps, nitmax, cdif):
    c = c0
    c_prev = c + cdif/2
    it = 1
    while not((np.abs(fun(c)) < eps and np.abs(c - c_prev) < cdif) or it > nitmax):
        c_prev = c
        c = c - fun(c) / dfun(c)
        it += 1
    return it, c, fun(c), c_prev, fun(c_prev), (c - c_prev)



In [16]:

    
newtraph2(2, f, df, my_stop, my_nitmax, my_cdif)









    Out[16]:





(7,
 3.1415926535897931,
 1.2246467991473532e-16,
 3.1415926536808043,
 -9.1011076112782062e-11,
 -9.1011198577461982e-11)



In [17]:

    
def secant2(a, b, fun, eps, nitmax, cdif):
    c = (a * fun(b) - b * fun(a)) / (fun(b) - fun(a))
    c_prev = c + cdif/2
    it = 1
    while not((np.abs(fun(c)) < eps and np.abs(c - c_prev) < cdif) or it > nitmax):
        a = b
        b = c
        c_prev = c
        c = (a * fun(b) - b * fun(a)) / (fun(b) - fun(a))
        it += 1
    return it, c, fun(c), c_prev, fun(c_prev), (c - c_prev)



In [18]:

    
secant2(2, 4, f, my_stop, my_nitmax, my_cdif)









    Out[18]:





(5,
 3.1415926535897931,
 1.2246467991473532e-16,
 3.1415926536048882,
 -1.5094913867333564e-11,
 -1.5095036332013478e-11)

Exercício 2

$f(x) = x + e^{-k x^{2}}$



In [107]:

    
from scipy.misc import derivative


def newtraphd(c0, fun, eps, nitmax):
    c = c0
    dfun = lambda x: derivative(fun, x, 0.0001)
    it = 1
    while np.abs(fun(c)) > eps and it < nitmax:
        c = c - (fun(c) / dfun(c))
        it += 1
    return it, c, fun(c), dfun(c)



In [97]:

    
f2 = lambda x, k: x + np.e ** (-k * x**2) * np.cos(x)



In [105]:

    
f2_k1 = lambda x: f2(x, 1)
newtraph(0, f2_k1, df2_k1, 1e-4, my_nitmax)









    Out[105]:





(15, -0.5883681212998364, 7.0172442686589065e-05)



In [122]:

    
f2_k50 = lambda x: f2(x, 50)
for i in range(1, 10+1): print(newtraphd(0, f2_k50, 1e-4, i))









    



(1, 0, 1.0, 0.99999999999988987)
(2, -1.0000000000001101, -1.0000000000001101, 0.99999999999988987)
(3, 1.1013412404281553e-13, 1.0000000000001101, 0.99999999998934275)
(4, -1.0000000000106573, -1.0000000000106573, 0.99999999999988987)
(5, 1.1013412404281553e-13, 1.0000000000001101, 0.99999999998934275)
(6, -1.0000000000106573, -1.0000000000106573, 0.99999999999988987)
(7, 1.1013412404281553e-13, 1.0000000000001101, 0.99999999998934275)
(8, -1.0000000000106573, -1.0000000000106573, 0.99999999999988987)
(9, 1.1013412404281553e-13, 1.0000000000001101, 0.99999999998934275)
(10, -1.0000000000106573, -1.0000000000106573, 0.99999999999988987)



In [123]:

    
for i in range(1, 10+1): print(newtraphd(-0.1, f2_k50, 1e-4, i))









    



(1, -0.1, 0.50350053278289919, 7.0955553139445682)
(2, -0.17095998981128269, 0.057578346024847515, 4.9465471855514309)
(3, -0.18260009832839749, 0.003047567616852459, 4.424209558854808)
(4, -0.18328893722183393, 1.0527388237396851e-05, 4.3936547460581012)
(4, -0.18328893722183393, 1.0527388237396851e-05, 4.3936547460581012)
(4, -0.18328893722183393, 1.0527388237396851e-05, 4.3936547460581012)
(4, -0.18328893722183393, 1.0527388237396851e-05, 4.3936547460581012)
(4, -0.18328893722183393, 1.0527388237396851e-05, 4.3936547460581012)
(4, -0.18328893722183393, 1.0527388237396851e-05, 4.3936547460581012)
(4, -0.18328893722183393, 1.0527388237396851e-05, 4.3936547460581012)